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We study the activation process in large assemblies of type II excitable units whose dynamics is influenced 
by two independent noise terms. The mean-field approach is applied to explicitly demonstrate that the assem¬ 
bly of excitable units can itself exhibit macroscopic excitable behavior. In order to facilitate the comparison 
between the excitable dynamics of a single unit and an assembly, we introduce three distinct formulations of 
the assembly activation event. Each formulation treats different aspects of the relevant phenomena, including 
the threshold-like behavior and the role of coherence of individual spikes. Statistical properties of the assembly 
activation process, such as the mean time-to-first pulse and the associated coefficient of variation, are found 
to be qualitatively analogous for all three formulations, as well as to resemble the results for a single unit. 
These analogies are shown to derive from the fact that global variables undergo a stochastic bifurcation from 
the stochastically stable fixed point to continuous oscillations. Local activation processes are analyzed in the 
light of the competition between the noise-led and the relaxation-driven dynamics. We also briefly report on a 
system-size anti-resonant effect displayed by the mean time-to-first pulse. 


I. INTRODUCTION 

For the wealth of local and intriguing collective phenomena 
displayed, the large assemblies comprised of excitable units 
have now been appreciated as a distinct class of dynamical 
systems. In terms of theory, the fundamental issue for under¬ 
standing these systems is whether their macroscopic dynamics 
itself exhibits excitable behavior. The other important issue 
naturally concerns the relation between such macroscopic ex¬ 
citability and noise. So far, stochastic effects have been iden¬ 
tified as a major factor contributing the collective behavior of 
systems of excitable units CHU. 

In this paper, we study large assemblies consisting of noisy 
type II excitable elements, which are represented by the 
canonical Fitzhugh-Nagumo ( FHN ) model. Conceptually, 
our focus will lie with two main points: (i) demonstrating 
that an assembly made up of excitable units can itself be con¬ 
sidered a macroscopic excitable element, and ( ii ) identifying 
the analogies and pointing out the differences between the ex¬ 
citable behaviors of a single unit and an assembly. 

Note that point (*), at least to our knowledge, has not been 
treated explicitly so far. In particular, the main obstacle for 
analytically approaching this issue is that the macroscopic dy¬ 
namics of an assembly of stochastic FHN units cannot be ex¬ 
pressed in a closed form via the global variables, which would 
otherwise make up a standard and the desired form of describ- 


* franovic@ipb.ac.rs 
I matjaz.perc@uni-mb.si 
i buric@ipb.ac.rs corresponding author 


ing the collective motion. Some alternative forms of analy¬ 
sis are not available due to complexity of the corresponding 
Fokker-Planck equation that assumes an integro-differential 
form. In order to resolve this and explicitly demonstrate the 
excitability feature at the assembly level, we use an approach 
that relies on the mean-field ( MF ) approximate model of the 
assembly’s collective motion. In several recent papers, the 
MF model has already proven successful in the analysis of 
systems described by large sets of stochastic (delay) differ¬ 
ential equations, in particular when treating stability of the 
stationary state, as well as the scenarios for the onset and the 
suppression of the collective mode lfl3l - fl5ll . 

Analysis on point (ii) can be carried out at two levels, one 
focused on the phenomenology involved, and the other con¬ 
cerning the statistical properties of the corresponding noise- 
driven activation processes. In terms of phenomenology, the 
excitability feature refers to capability of systems to generate 
spiking (pulse-like) responses or small-amplitude excitations, 
which are separated by some form of threshold. For a single 
unit, the large-amplitude response is comprised of the activa¬ 
tion and the relaxation stage, such that the former is strongly 
influenced by noise, whereas the latter is typically determinis¬ 
tic and maintains a stereotype profile in a broad range of noise 
values. Among else, the stereotype character of pulse implies 
that its amplitude and width are independent on the form of 
perturbation applied. 

The stated arguments on the notion of excitability should 
naturally hold in case of an assembly as well. Nevertheless, 
what may be different are the details related to the assembly’s 
threshold-like behavior, which by itself stands out as a highly 
non-trivial issue. Another point of difference concerns the lo¬ 
cal mechanisms by which excitations of units within a popula- 
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tion are elicited. In particular, each unit may be evoked to emit 
a pulse either by noise or via the interaction terms. Adhering 
to formulation of activation event for an excitable unit stated 
in our previous paper m, the latter point does not merely 
imply that the activation processes of individual units cannot 
be considered independent, as if they were driven by uncor¬ 
related random perturbations, but more significantly indicates 
that the activation process of an arbitrary unit may be influ¬ 
enced by the relaxation processes of other units. The third 
point one should stress concerns the process of spike emis¬ 
sion at the assembly level. In particular, having accepted the 
description of collective motion in terms of global variables, 
one should also bear in mind that their amplitude is affected 
by coherence of individual spikes, such that the increased co¬ 
herence implies larger amplitude of the relevant macroscopic 
variable. By such a concept, evoking a large-amplitude exci¬ 
tation requires that the spikes for a sufficient fraction of indi¬ 
vidual units become coherent. Analyzing coaction of activa¬ 
tion and relaxation processes for individual elements, as well 
as the role of spike coherence with respect to assembly pulse 
emission, present intricate issues, absent if small groups of 
units are considered. 

Apart from phenomenological aspects, the comparative 
analysis of excitable behaviors of a unit and an assembly will 
address the respective noise-driven activation processes. In 
analogyto our paper on a single and two interacting excitable 
units j21J], we associate the term activation solely to the as¬ 
sembly’s large-amplitude excitation. In principle, the defini¬ 
tion of assembly activation problem should incorporate appro¬ 
priate boundary conditions, which provide a clear-cut distinc¬ 
tion between the small- and the large-amplitude excitations. 
Nevertheless, while a single most adequate and relevant defi¬ 
nition for the activation problem may be given in cases of an 
excitable unit or a pair of units |21], the specific character of 
global variables, or rather the fact that their amplitudes depend 
on coherence of individual spikes, prevents us from establish¬ 
ing such a formulation in case of an assembly of excitable ele¬ 
ments. Instead, we introduce three alternative formulations of 
the assembly activation event which emphasize different as¬ 
pects of the relevant phenomena. Their merit will depend on 
the aims of the particular study of assembly dynamics and the 
fashion in which one can adapt them to potential applications. 

Two of the formulations rest on the standard description 
of collective motion in terms of global variables, and are in¬ 
tended precisely at examining the analogy between the ex¬ 
citable behaviors of a single unit and an assembly. In par¬ 
ticular, one formulation is consistent the threshold boundary 
approach for a FHN unit (terminating boundary condition 
given by a threshold x-value), whereas the other derives from 
characteristic boundary approach (terminating boundary con¬ 
ditions given by an appropriate boundary set). Nevertheless, 
for comparison we also adopt a formulation where the assem¬ 
bly activation is treated as a compound event, comprised of 
only first-pulse responses of a sufficient fraction of participant 
units, regardless of whether the emitted spikes are coherent or 
not. The implications of the three formulations are analyzed in 
terms of statistical properties of the corresponding activation 
events, characterized by the dependencies of the time-to-first 


pulse emission ( TFP ) and the related coefficient of variation 
on noise. 

The paper is organized as follows. In Sec. [TT] is provided 
the background on the applied model. Section Hill addresses 
the details of the assembly’s excitable behavior, including the 
analysis carried out on the deterministic MF approximation. 
In Section||V]are laid out the three formulations of the assem¬ 
bly activation event. Section [V| contains the detailed numer¬ 
ical analysis on the statistical properties of activation events 
conforming to the three adopted formulations. We also con¬ 
sider the qualitative explanation for the bimodal or unimodal 
distributions of local activation events typical for certain do¬ 
mains of noise intensities. Subsection IV Dl concerns the ef¬ 
fects related to system size, including the ’’anti-resonance” 
found for the mean TFPs at fixed noise intensities. In Sec. 
IVIl is given a brief summary of our main points. 


II. DETAILS OF THE APPLIED MODEL 

As a paradigm for analyzing collective excitable behavior, 
we consider an assembly comprised of FHN units. The dy¬ 
namics of an arbitrary unit i is given by 

N 

dxi = (Xi — rcf/3 — yi)dt + ^(xj; - Xi)dt + y/2DidW[ 

JV i =i 

dyi = e(xi + b)dt + \j2DidW\, i = 1,.. .N (1) 

where b and e are the intrinsic unit parameters, while c denotes 
the coupling strength. The units are assumed to be identical 
and are connected in the all-to-all fashion. Parameter e is set 
to a small value e = 0.05, which warrants that the character¬ 
istic time scales for Xi and y, evolution are sharply separated. 
Being type II excitable means that the units are poised close to 
transition toward oscillatory state via the Hopf bifurcation |2ll . 
The bifurcation parameter b is set to 1.05, the value just below 
critical threshold. Excitability feature of a single FHN unit 
has been extensively analyzed lli 12}, 120], and an overview can 
also be found in our preceding paper |21]. At variance with 
the latter, the perturbation here may either arrive from the in¬ 
teraction terms, or may be caused by random fluctuations due 
to two independent sources of noise. Motivated by the pos¬ 
sible interpretation in the field of neuroscience fl9l l20h . we 
adopt the convention by which the stochastic terms in the fast 
(slow) variables are referred to as external noise D\ (internal 
noise Z) 2 ). Note that dW£, k £ {1, 2}, i = 1 ,N denote 
stochastic increments of independent Wiener processes whose 
averages and correlations satisfy (dWJ.) = 0, (dW'{.dWj l ) = 
dtSkiSij. 

In order to gain insight into the assembly’s collective dy¬ 
namics, one may first carry out the bifurcation analysis of 
the deterministic (noiseless) version of system ([]}• For N 
sufficiently large so that the terms 0((c/N) 2 ) and of higher 
order can be neglected, the characteristic equation describ¬ 
ing the stability of equilibrium (xi,yi, x 2 , t/ 2 , • ■ •, Xjv, t/iv) = 
(—5, —b + b 3 / 3, —6, —b + b 3 / 3,..., —b, —b + b 3 / 3) is given 
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by an approximate expression 

(A 2 - (1 - b 2 ) A + e)(A 2 - (1 - b 2 - c)A + e)^- 1 = 0. (2) 

Since b = 1.05 is kept fixed, it follows that the equilibrium 
may become unstable only via the direct supercritical Hopf 
bifurcation controlled by c. Nonetheless, the critical c value 
is c* = 1 - b 2 < 0, which implies that the positive cou¬ 
pling strengths do not affect the stability of equilibrium. In the 
present paper, we only consider the subcritical values c > 0, 
such that the system 0 always lies in the excitable regime. 

As for the impact of stochastic fluctuations on the asymp¬ 
totic collective dynamics, it is known that the noise intensity 
may act as a control parameter in excitable media, giving rise 
to three generic regimes of macroscopic behavior y|]. In par¬ 
ticular, when noise is systematically increased, the dynam¬ 
ics of global variables undergoes a sequence of transitions, 
first exhibiting the stochastically stable equilibrium, then the 
stochastically stable limit cycle, and eventually it decays into 
disordered behavior. This sequence may be explained as fol¬ 
lows. In the approximately stationary state, at any given mo¬ 
ment, most of the units lie in the vicinity of equilibrium, 
whereas the relatively rare excursions due to weak noise or the 
interaction terms remain incoherent. Therefore, the macro¬ 
scopic variables are only marginally different from equilib¬ 
rium values of single elements. At some point, the increase 
of noise induces more or less coherent oscillations of units 
which are easy to synchronize. This conforms to the onset 
of the collective mode according to the scenario of stochas¬ 
tic bifurcation li22l - t25ll . In the supercritical state, the global 
variables follow a limit cycle attractor whose profile is simi¬ 
lar to that of relaxation oscillations of individual units. Once 
the noise intensity becomes strong enough to overcome the 
effect of couplings, the disordered regime sets in. The spik¬ 
ing frequency of units remains high, but their activity desyn¬ 
chronizes. Since the majority of units at any moment are re¬ 
fractory, global variables display irregular oscillations with a 
quenched amplitude. 

From our perspective, the conceptually most important 
transition is the one associated to occurrence of noise induced 
collective oscillations. This phenomenon is relevant for the 
activation process because it indicates the loss of stochastic 
stability for the fixed point. In other words, under increas¬ 
ing noise, the attractive ability of the fixed point gradually 
reduces and is eventually lost, which naturally affects how the 
phase point escapes the vicinity of equilibrium. Note that the 
above sequence of states has previously been verified in case 
of an assembly of FHN elements driven by internal noise. 
We have found that the similar sequence persists if the units 
are influenced by external noise, though the individual activ¬ 
ities become more difficult to synchronize, rendering the am¬ 
plitude of collective oscillations comparably smaller than the 
one emerging under internal noise. 


III. EXCITABLE BEHAVIOR AT THE ASSEMBLY LEVEL 

Having summarized the points relevant for the determin¬ 
istic part of dynamics described by 0, we turn to charac¬ 


terization of the assembly’s excitability feature. The stan¬ 
dard approach to collective motion is to introduce the macro¬ 


N 


scopic variables X(t) = ( Xi(t )) = j? %i(t) and Y(t) = 


i-1 


N 


(yi(t )) = X ^2 yiit), whereby the aim typically lies in es- 

2=1 

tablishing some form of analogy between the dynamics of 
single units and an assembly. The latter would be espe¬ 
cially relevant for examining the issue of excitability at the 
assembly level. Nevertheless, it is evident that the com¬ 
pound effect of nonlinear terms prevents the whole assem¬ 
bly to evolve in the fashion analogous to that of a single unit, 
which would occur only if (x 3 ) = {xi) 3 were to hold. There¬ 
fore, the collective dynamics in principle cannot be expressed 
in a closed form via the global variables. The other poten¬ 
tial approaches to analysis of macroscopic excitable behavior 
are severely limited by the difficulties associated to Fokker- 
Plank formalism, where the equation for the one-particle den¬ 
sity P(x, y, t) acquires a complex integro-differential form 
■§tP = --jk[ x (X- c )-^--y+ c fxiP(x 1 ,y 1 ,t)dx 1 dy 1 }P- 

o. e{x + b)P + Dl ^P +D 2 ^P. 

The arguments above imply that one is required to intro¬ 
duce approximations for collective dynamics corresponding 
to the deterministic part of system 0 in order to analyze 
the assembly excitability feature and the associated threshold¬ 
like behavior. In the following, we consider two approxi¬ 
mate models of collective motion, distinguished by the fash¬ 
ion in which the effect of interaction terms is resolved. The 
first model holds if the interaction terms vanish or can be ne¬ 
glected. Note that this condition is satisfied if the initial con¬ 
ditions for all the units are identical or lie close to each other. 
Since the evolution of the system is deterministic and the cou¬ 
pling is diffusive, such selection of initial conditions facilitates 
that the whole assembly acts as a macroscopic excitable ele¬ 
ment. The evolution of global variables is then given by the 
equations analogous to those for a single unit 


dX = {X- X 3 /3 - Y)dt 
dY = e(X + b)dt. 


(3) 


The details regarding the excitability feature of such a system 
are well established |20j. In particular, recall that its threshold 
behavior is associated to ’’ghost-separatrix”, a thin layer made 
up of canard-like trajectories that foliate around the maximum 
of the fast-variable nullcline, whereby the spread increases 
with the characteristic scale separation ratio e. 

Let us now consider a more sophisticated approximation 
that takes into account the net effect of the interaction terms. 
The analytical framework suitable for demonstrating the as¬ 
sembly excitability feature in this more general case is pro¬ 
vided by the mean-field approach. Before laying out the de¬ 
tails, we make an overview of the ingredients crucial for the 
derivation of the MF model, as well as the results achieved so 
far on treating the systems of large sets of stochastic fdelay)- 
differential equations. In principle, deriving the MF model 
involves a number of nontrivial elements, and it ultimately 
leads to a deterministic system amenable to standard bifur¬ 
cation analysis, where noise intensity may act as a bifurcation 
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parameter. The MF method combines the cumulant approach 
with Gaussian approximation, according to which all the cu- 
mulants above second order are assumed to vanish. The latter 
is intended as a closure hypothesis, which is necessary due 
to presence of nonlinear terms in the exact system 0. Thus, 
starting from the original system which in general comprises 
kN (delay) differential equations, where k is the number of 
local degrees of freedom, one ends up with a set of k(k + 3)/2 
deterministic (delayed) equations describing the evolution of 
the means, as well as the appropriate variances and the covari¬ 
ances. 

As for the main results achieved thus far, the MF method 
has already been applied in analyzing the stability of assem¬ 
blies of (delay)-coupled excitable elements, as well as the 
scenarios for the onset and the suppression of the collective 
mode. In particular, the bifurcations displayed by the MF 
model can qualitatively account for the stochastic bifurcations 
which the exact system undergoes H1 i S El EMl- 
it has also been shown that the MF model can provide ac¬ 
curate quantitative predictions, reflected in a close agreement 
between the oscillation period of the MF model and the av¬ 
erage oscillation period of the exact system El El- Further¬ 
more, the MF approach has proven successful in the analysis 
on stability and the onset/suppression of the collective mode 
in case of interacting assemblies, thereby indicating a poten¬ 
tial extension to modular networks El- We stress that the 
approximations behind the MF model used here, called the 
Gaussian approximation and the quasi-independence approx¬ 
imation, have been analyzed in detail, not only in terms of 
precise formulations and adaptation for the systems of class 
II excitable units, but also with respect to parameter domains 
that warrant their validity El . In this context, an important 
finding is that the dynamics of the MF model can indicate 
in a self-consistent fashion the parameter domains where the 
MF approximations break down. 

The MF model corresponding to 0 involves five equa¬ 
tions describing the evolution of the means m x (t ) = 

E(xi(t)),m y (t) = E(yi(t)), the variances s x (t) = 

E((xi(t) - ( Xi(t))) 2 ),s y (t ) = E((yi(t ) - (t/i(f))) 2 ) and the 
covariance u(t) = E((xi(t)-(x i (t)))(y i {t)-(y i (t)))). Note 
that E(-) denotes expectations over an ensemble of different 
stochastic realizations. The detailed derivation of the MF 
model may be found in El El . whereas here we just state 
the final result 

ni x = m x m x (t ) 3 /3 - s x m x - m y , 

m y = e(m x + b), 

s x = 2s x (l — m x — s x — c) — 2u + 2D\ 
s y = 2 eu + 2 D 2 , 

u = u(l - ml - s x — c) + es x - s y . (4) 

As already indicated, the influence of stochastic terms is ex¬ 
pressed through the noise intensities 1) \ and D^. 

Nevertheless, in order to make the analogies between the 
excitable behaviors of a single FHN unit and an assembly 
explicit, one should arrive at the system describing the collec¬ 
tive motion by two equations. The latter would allow one to 
apply the phase plane analysis derived from the framework of 


singular perturbation theory. This is achieved by introducing 
an additional ’’adiabatic” approximation Eh, which consists 
in assuming that the relaxation of second-order moments is 
much faster than that of the first-order ones. This is not a 
crude approximation given that the initial conditions of units 
in the exact system are set to be identical (coinciding with 
the deterministic fixed point), especially if the noise intensi¬ 
ties are not too large. Having replaced the fast variables with 
the stationary values, one obtains the following system for the 
dynamics of the means: 


dm x 

dt 


= F(m x ,m y ) =m x - -m x - m y -—(1 

ml + y/(l — c— ml) 2 + 4(Di + D 2 /e) 


— c— 


dm y 

dt 


G(m x ,m y ) = e(m x + b). 


(5) 


An apparent advantage of ([5]) is that one may extend the re¬ 
sults of phase plane analysis to assembly dynamics, thereby 
gaining insight into whether and how the excitability feature 
is manifested at the level of global variables. 

The main point is that the system (0 displays class II ex¬ 
citable behavior, which provides an indication that the collec¬ 
tion of excitable FHN units described by 0 itself constitutes 
a macroscopic excitable system. The m x and m y nullclines, 
as well as an illustration of how the two marginally differ¬ 
ent initial conditions give rise to small- or large-amplitude re¬ 
sponses are provided in Fig. E] Note that the results are ob¬ 
tained for system 0 under Di = D 2 = 0. The m, T -nu lie line 
again consists of three branches, such that in the singular limit 
e — > 0 the spiking branch Ss and the refractory branch Sr are 
attractive, whereas the middle branch is unstable. Compared 
to the case of single unit, the profile of the middle branch is 
changed and includes a flexion point, which may be attributed 
to the compound effect of interaction between the units. The 
two types of population response to perturbation, or rather the 
associated threshold-like behavior, imply the existence of a 
soft boundary between the corresponding initial conditions. 

Carrying out the analysis analogous to that for a single unit, 
described in brief in the caption of Fig. E] one may show that 
the boundary is again given by a thin layer of system trajecto¬ 
ries which foliate around the maximum of the m.,.-nuHcline. 
The relevant segments of such trajectories are indicated in Fig. 
Q]by the black solid lines. What is found is quite reminiscent 
to ’’ghost-separatrix” in case of a single excitable unit (20]. 
The foliation here also becomes more pronounced with in¬ 
creasing e. Compared to those for a single unit [20], the tra¬ 
jectories that make up the boundary layer are seen to converge 
further away from the maximum, mainly because the unstable 
branch in case of an assembly involves a flexion point. Also, 
the numerical evidence suggests that the mean-field variables, 
and therefore the collective dynamics of the assembly, shows 
greater sensitivity to initial conditions compared to that of a 
single unit. This conclusion follows from the fact that under 
analogous parameter values, obtaining the trajectories that be¬ 
long to boundary set requires a higher numerical precision in 
case of the MF model than for a single unit, cf. [21] and Fig. 

El 
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FIG. 1. (Color online) Characterization of the assembly excitable be¬ 
havior via the phase plane analysis of the corresponding MF model 
The equilibrium (EQ) lies at the intersection of the nullclines 
F(m x ,m y ) = 0 and G(m x ) = 0. The mz-nullcline is comprised 
of three branches, whereby the spiking and the refractory branches 
Ss and Sr are attractive. For hnite e, the boundary between the 
sets of initial conditions that lead to small- or large-amplitude re¬ 
sponses, SAR and LAR respectively, foliates into a thin layer of 
canard-like trajectories ( CNR ), which we refer to as ’’ghost separa- 
trix” (solid black lines). The trajectories belonging to the boundary 
layer are obtained by fixing the m, initial condition to a particular 
value m x fi = —3, while sweeping over m y fi. The fact that a differ¬ 
ence in rriyfl of the order of 10 -18 evokes a different type of response 
corroborates extreme sensitivity to initial conditions in vicinity of 
ghost separatrix. The results shown refer to case t = 0.07, b = 1.05. 


Note that in Sec. 0 we use the MF model to gain insight 
into the stochastic stability of equilibrium of the exact sys¬ 
tem, or rather the latter’s sensitivity to random perturbations. 
This point will be crucial for explaining how the form of the 
t{D\, D 2 ) dependence is related to the stochastic bifurcation 
leading from the stochastically stable fixed point to continu¬ 
ous oscillations. 


IV. ALTERNATIVE FORMULATIONS OF ASSEMBLY 
ACTIVATION PROBLEM 

As announced in the Introduction, we consider three alter¬ 
native formulations of the activation problem, which are asso¬ 
ciated to different definitions of the assembly activation event. 
Since the conceptual aspects of a large assembly’s activation 
problem have not been treated so far, one cannot a priori hold 
any formulation preferred over the others. Thus, the impli¬ 
cations of each of the formulations will be qualitatively ana¬ 
lyzed and then compared. Note that the physical picture laid 
out here is distinct from the cases of a single and two cou¬ 
pled excitable units treated in our previous paper [21], where 
we have been able to provide a unique adequate formulation 
of the activation event, which among else, allows an immedi¬ 


ate generalization from a single unit to a two-unit setup. To 
properly address the issue of assembly activation, one should 
invoke a couple of remarks from that study. 

First, the activation problem we consider is associated to 
pulse emission, and as such cannot be viewed as an extension 
of a typical escape problem. Note that the latter would require 
a genuine saddle structure at the terminating boundary |[26], 
which implies coexistence between two attractors, whereas in 
our problem the fixed point is the only relevant, and often the 
unique attractor. In terms of specifying the terminating bound¬ 
ary, an important point has been to replace the somewhat arbi¬ 
trary notion of ’’threshold” for pulse emission by the relevant 
boundary set consistent with the underlying structure of the 
phase space. Finally, in case of two units, we have argued 
that the definition of activation event where the phase points 
of each unit are supposed to reach the appropriate terminat¬ 
ing boundary is preferred over any formulation involving the 
two-unit averages (xi(t) + X 2 {t ))/2 and (yi(t) + 2 / 2 (*))/2. 
This applies because the averages contain additional informa¬ 
tion on synchronization of units, which are secondary to the 
two-unit activation process. Nevertheless, such an approach 
can only be maintained for small groups of units. For larger 
assemblies, it is of interest to state the activation problem in 
terms of global variables, since they present a standard tool 
for describing collective motion, both theoretically and in ap¬ 
plications. 

The point which makes the case of large populations in¬ 
triguing concerns the underlying mechanisms of activation. In 
particular, the processes on a local and global level are not in¬ 
fluenced only by noise, but are also strongly contributed by 
the relaxation of units. The arguments above further suggest 
that an elaborate study on analogies and differences between 
the statistical properties of activation process in small groups 
of units and the large assemblies would be in order. Therefore, 
our approach is on one hand to retain the formulation of the 
assembly activation problem inherited from the two-unit case 
where the global variables are not considered, and on the other 
hand, to introduce two additional formulations that explicitly 
refer to dynamics of global variables. These three formula¬ 
tions are specified as follows. 

Formulation 1 .The assembly activation event occurs when 
more than a half of participant units have emitted their first 
pulses. According to this, assembly activation is perceived as 
a compound event made up of local activation events. For the 
local events we adopt the definition provided in our previous 
paper |21]: the activation path of a single excitable unit in¬ 
fluenced by noise emanates from the deterministic fixed point 
and terminates at the boundary set coinciding with the spik¬ 
ing branch of the limit cycle which would exist in the corre¬ 
sponding supercritical state. In the present context, the termi¬ 
nating boundary set can with sufficient accuracy be approxi¬ 
mated by the spiking branch of the cubic nullcline for a non¬ 
interacting unit. Motivation behind formulation 1 draws in 
part from certain applications, especially in the field of neu¬ 
roscience, where population response to external stimuli typi¬ 
cally engages a certain fraction of units, rather than the entire 
assembly Il27l432ll . From the qualitative perspective, demand¬ 
ing any reasonable macroscopic fraction other than a half of 
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units to be activated makes as good a choice as any, because 
the main statistical properties of the ensuing activation pro¬ 
cess will remain similar. In a sense, formulation 1 can be 
interpreted as an extension of the definition introduced for a 
two-unit activation event in our previous paper (21]. 

Formulation 2.The assembly activation event occurs if 
the global variable X(t) crosses the predefined threshold 
Xp(X(t) > Xq. X'(t) > 0). Unlike the case of a single 
unit I20l I2ltl . the formulation involving an explicit threshold 
is justified for the global variable because its amplitude de¬ 
pends on coherence of spikes of single units. The latter point 
introduces ambiguities when attempting to analyze the assem¬ 
bly threshold-like behavior. In other words, for a single unit, it 
is not difficult to distinguish between the small and large am¬ 
plitude responses, given that the amplitude of superthreshold 
excitations is stereotypical. However, in case of an assembly, 
one is able to understand the associated threshold-like behav¬ 
ior only in terms of the approximate MF model, cf. Sec. |III1 
but cannot provide clear-cut criteria in terms of global vari¬ 
ables of the exact system, especially if coherence of units’ ac¬ 
tivities for the given parameter set is weak. This point will be 
further explained when discussing formulation 3, while here 
we just mention that noise domains may be found where for¬ 
mulation 2 is more or less suitable. Compared to formulation 
1, formulation 2 is conceptually distinct because the assembly 
activation events can be affected by multiple spikes of individ¬ 
ual units. The potential consequences of selecting particular 
threshold values Xq will be discussed in the next Section. 

Formulation 3.The assembly activation event occurs once 
the phase point associated to global variables reaches the ap¬ 
propriate terminating boundary set, defined in analogy to a 
single unit case. This formulation involves an implicit as¬ 
sumption that the threshold-like behavior of an assembly is 
qualitatively similar to that of individual units. As indicated 
in Sec. mn if the condition (x 3 (t)) ~ ( Xi(t )} 3 is fulfilled, 
the dynamics of global variables can be approximated by the 
equations analogous to that of a single unit. Consistent with 
this approximation, the collective motion of finite assemblies 
in presence of noise may be described by ll3[|3jl 


dX = (X - X 3 /3 - Y)dt + yj‘^-dW 1 
/ o n 

dY = e(X + b)dt + y -^-dW 2 . (6) 

As a corollary, one may use the phase plane analysis and de¬ 
scribe the assembly first-pulse emission in terms of motion of 
phase point ( X , Y) to the spiking branch of the corresponding 
A’-nu lie line, which is then naturally considered the relevant 
terminating boundary set for the assembly activation problem. 
Nevertheless, the given physical picture is valid only if the 
fluctuations of individual units are small enough for the above 
condition to apply, which ultimately depends on the noise in¬ 
tensities. 

An important note on formulation 3 is that coherence be¬ 
tween spikes (approximate matching of spike times) of a large 
fraction of individual units is required to elicit an assembly ac¬ 
tivation event, viz. the first-pulse emission for the global vari- 



X 

FIG. 2. (Color online) Illustration of differences between formula¬ 
tions 2 and 3. For certain ( Di,D 2 ), the assembly threshold-like 
behavior is difficult to resolve because the amplitude of large fluc¬ 
tuations is not stereotypical. The presented (X, Y ) orbit is ob¬ 
tained for a single stochastic realization under parameter set D\ = 
0.0004, D 2 = 0.0008, c = 0.1, N = 50. The bold dotted line in¬ 
dicates the .Y-nullcline consistent with the approximate system 
Apart from small-amplitude excitations that remain in close vicinity 
of deterministic fixed point, one also finds substantially larger exci¬ 
tations far above the stochastic fluctuations around the "stable state”, 
where the phase point may either fall short of the spiking branch of 
the A'-nullcline (segment of trajectory shown by the orange dashed 
line), or may actually reach it (segment indicated by the solid black 
line). The former instance conforms to activation event by formu¬ 
lation 2, but does not comply with formulation 3. The criteria as 
to which large excitation is considered an assembly activation ulti¬ 
mately depends on the scope of the particular study. 


ables. For this point, let us make additional remarks on the is¬ 
sues of terminating boundary conditions and the role of spike 
coherence in the assembly activation process. First, there is no 
principle difference in the dynamics underlying activation sce¬ 
narios for three formulations of the assembly activation event. 
The differences in statistical features, which will be discussed 
in the next Section, are caused by the selection of terminating 
boundary conditions. This is made explicit in Fig. [2] where 
a single stochastic orbit for ( X, Y) is used to illustrate the 
difference between formulations 2 and 3. According to for¬ 
mulation 2, the assembly has emitted a pulse once X crosses 
a certain threshold (segment of the trajectory indicated by the 
dashed red line), while formulation 3 requires that the phase 
point reaches the spiking branch of the X-nullcline (segment 
of the trajectory indicated by the solid black line). Thus, con¬ 
sidering the very same stochastic trajectory for (X, Y), the 
assembly activation event conforming to formulation 2 would 
have happened earlier than the one satisfying formulation 3. 
Nevertheless, both of the considered sections of the (A', Y) 
orbit correspond to large excitations far above the stochastic 
fluctuations in vicinity of equilibrium. In fact, which excita¬ 
tion should be considered an activation event effectively de¬ 
pends on the scope of the study. In this context, formulation 2 
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FIG. 3. (Color online) Mean TFP s for the assembly activation process and the role of stochastic bifurcation, (a), (b) and (c) show r(D i, D 2 ) 
dependencies for formulations 1 — 3 of the activation event, respectively. Results are obtained for c = 0.1 and the assembly size N = 100, 
whereby the stochastic average is taken over an ensemble of 300 different activation paths. The inset in (a) illustrates exponential divergence 
of r(Z? 2 ) for fixed D\ = 10 -5 as one approaches D% « 1.25 * 10 -5 (the value indicated by the dotted line). The qualitative similarity 
between the t(D 1 , D 2 ) dependencies in (a), (b) and (c) is associated to stochastic bifurcation underlying transition to continuous oscillations. 
The corresponding bifurcation curve D 2 (Di), determined by analyzing MF model (5), is shown in (d). 


is more adapted to practical applications, whereas formulation 
3 has a more theoretical background. 

Comparing Fig. Q]and Fig. [2j one should point out that the 
latter captures certain fine details on collective dynamics that 
cannot be reflected in the former. This is because the physical 
background of Fig. |T] lies in the approximate MF model, 
where the stochastic fluctuation effects above the second order 
are averaged out. In other words, the behavior of the exact 
system in Fig. [2] is associated with large fluctuations, which 
can no longer be described by the MF model. 


V. STATISTICAL PROPERTIES OF ACTIVATION 
PROCESS 

A. Mean TFPs and TFP variability 

The statistical properties of the assembly activation pro¬ 
cess influenced by external and internal local noise are char¬ 
acterized by the mean TFP t(Di , D?) and the associated 
coefficient of variation R(Di,D 2 ) Il34l436il . Both quanti¬ 
ties involve averaging over an ensemble of different stochas¬ 
tic realizations. In particular, the mean TFP is given by 


lL r 

t(D 1j D 2 ) = ( T k ) = 7 T £ T fe (Di,D 2 ), where n r denotes 
r k =1 

the number of realizations, while the expression for the coef¬ 
ficient of variation reads R{D\,D 2 ) = '^ < ' Tk { Tk ) T ~- Since 
R is variation of TFPs normalized by the mean, its smaller 
values indicate a better clustering of individual TFP s around 
t m. Note that the numerical simulations are carried out via 
Heun integration scheme with the fixed time step St = 0.002. 
The results for the mean TFPs and the variances are ob¬ 
tained by averaging over an ensemble of at least 300 different 
stochastic realizations of the activation process. For each re¬ 
alization, the initial conditions for all units are identical and 
are given by the deterministic fixed point of system ©■ 

The fields t(D±, D 2 ) corresponding to formulations 1 — 3 
of the assembly activation event are plotted in Figs. [3J a ), EJb) 
and |3jc), respectively. Note that all three dependencies ex¬ 
hibit qualitatively similar behavior, whereby one can clearly 
discern the domain of large mean TFPs and the plateau re¬ 
gion. We find that these analogies derive from the fact that the 
profile of t(D\, D 2 ) dependencies in each of the instances is 
crucially influenced by the stochastic bifurcation. The latter 
corresponds to the noise-induced transition from the stochas¬ 
tically stable fixed point to continuous oscillations. It is in- 
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FIG. 4. (Color online) Build-up of spike coherence prior to assembly 
activation event. The data concern a single stochastic realization of 
the assembly activation event consistent with formulation 2 („Yo = 
0.4). The dotted line shows the X(t) dependence, whereas the solid 
lines refer to the appropriate time series of individual units Xi(t). 
The system parameters are D\ = 0.0001365, D 2 = 0.0001, c = 
0.1, IV = 100. 


tuitively clear that the stochastic bifurcation should be asso¬ 
ciated to boundary between the two characteristic forms of r 
behavior, because above the bifurcation, the fixed point loses 
its attractive power, which makes it easier for noise to in¬ 
duce large-amplitude fluctuations. Also, once the attractive 
power of the fixed point is lost, any further increase of noise 
cannot induce qualitatively novel effects. This fact accounts 
for the existence of the large plateau region in the t(D\, I)2) 
dependence. As already indicated, we may gain insight into 
stochastic bifurcation by carrying out bifurcation analysis on 
the MF counterpart (0 of the exact system. In particular, it is 
demonstrated that the MF model undergoes direct supercriti¬ 
cal Hopf bifurcation which qualitatively reflects the stochastic 
bifurcation of the exact system. The Hopf bifurcation curve 
D 2 {D\) obtained for the MF model at fixed c = 0.1 is shown 
in Fig. 0d). 

Note that for small (1) \. D 2 ) the mean TFP s grow ex¬ 
tremely large. This is explicitly illustrated in the inset of Fig. 
0 a), which shows shows that r(£> 2 ) (for fixed small 1) \) ex¬ 
ponentially diverges as a certain value of noise intensity I)* 2 
is approached. The latter point clearly demonstrates the ex¬ 
istence of potential barrier associated to the assembly’s acti¬ 
vation process. The lower boundary on noise intensities that 
facilitate activation is the largest for formulation 3, and de¬ 
creases for formulations 2 and 1. The reason behind such a 
behavior will be considered below. 

In quantitative terms, the r dependencies for formulations 
1 — 3 are manifestly different in the regime subcritical rela¬ 
tive to stochastic bifurcation. In particular. Fig. 0 shows that 
the corresponding r values for the same noise intensities are 
typically ordered as t\ < T 2 < 73 . This maybe explained by 
looking into the role played by the spike coherence between 



FIG. 5. Impact of different Xo on the results for formulation 2 of ac¬ 
tivation event. The example provided illustrates the form of t{Xq) 
dependence typical in a broad range of (D\. D 2 ,c : N) values. One 
finds similar behavior for all the other statistical quantities consid¬ 
ered. The particular data are obtained for Di — 0.0207018, O 2 = 
0.00101739, c = 0.1, N = 100, with the stochastic averaging car¬ 
ried out over an ensemble of 500 different activation paths. The re¬ 
sults imply that the statistical properties of activation process derived 
from formulation 2 are qualitatively independent on Xq. 


the individual units in the respective assembly activation pro¬ 
cesses. Let us focus first on the adjustment between the local 
time series Xi(t) prior to assembly activation, see Fig. 0 Nat¬ 
urally, the early spikes of single units will be the least coherent 
because the firing is driven by noise, and the effect of diffusive 
coupling in bringing their firing times closer is weakly felt. 
As the time passes, the impact of coupling on units’ dynamics 
is felt more strongly, which gradually leads to a build-up of 
spike coherence within the assembly. In other words, firing of 
individual units will become more coherent (more spikes fall 
within a relatively narrow time window) as the result of in¬ 
teractions. With time, such approximate synchronization may 
comprise an increasing number of units, depending on the c 
value. 

Now we consider how this reflects to assembly activation 
for different formulations of the activation event. Formulation 

1 requires only that more than half of units have emitted their 
first spikes, imposing no demand for their coherence. In line 
with the above arguments, satisfying such a requirement will 
take the least amount of time, which means that the pertaining 
r will be smaller compared to mean TFPs for formulations 

2 and 3 under the same parameter set. Formulations 2 and 3 
concern the global variables and thereby do require a certain 
level of spike coherence within the assembly. Formulation 3 
imposes a more strict condition, so its fulfillment is associ¬ 
ated to the largest level of spike coherence. Achieving larger 
level of spike coherence takes more time, so in the subcriti¬ 
cal regime TFP s for formulation 3 will necessarily be larger 
than those for formulation 2. It is also apparent that the assem- 
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FIG. 6. Impact of two noise sources on variability of assembly activation process. The R values corresponding to formulations 1 — 3 of the 
activation event are represented by open squares, solid circles and open triangles, respectively. In (a) is provided the R(D i) dependence for 
fixed Z >2 = 0.00005, whereas (b) shows the R(D 2 ) behavior at fixed D\ = 0.00062. (c) refers to R(Di ) dependence for D 2 = 0.00459. 
The noise values kept fixed in (a) and (b) are small, such that the system undergoes stochastic bifurcation under variation of D 1 in (a) and D 2 
in (b). The D 2 value in (c) is supercritical relative to stochastic bifurcation. The remaining system parameters are c = 0.1, N = 100. 


bly activation will necessarily be induced by the first spikes of 
single units only for formulation 1, whereas activation by for¬ 
mulations 2 and 3 is bound to be contributed by latter spikes of 
individual units. Above the stochastic bifurcation, coherence 
of spikes between the units is naturally achieved, such that the 
assembly activation for all three formulations is contributed 
by the first spikes of units. This is the reason why the mean 
TFPs for (D 1 , D 2 ) above the bifurcation have similar values 
for all three formulations of the activation event. 

Before proceeding, we make a brief remark on certain de¬ 
tails related to formulation 2 of the activation event. As indi¬ 
cated in Sec. [IV] it has to be verified whether the results on 
statistical properties of the activation process depend on the 
particular value of the threshold A'o- We have established that 
for any reasonable selection, viz. the Xo values lying suffi¬ 
ciently above the amplitude of stochastic fluctuations around 
the fixed point, the corresponding dependencies t(1)\ . I)->) 
maintain qualitatively similar profile. This holds for X (J val¬ 
ues within the range Xj £ [—0.1,1.3]. Considering r val¬ 
ues as an example, we have shown that r(Xo) for arbitrary 
, D 2 ) monotonously increases until reaching saturation 
around Xo « 1.2, see Fig. 0 The analogous conclusions 
regarding the qualitative independence of the results on Xo 
have been confirmed for all the other statistical properties of 
the activation process. 

Let us now examine the impact of two different noise 
sources on the behavior of R in view of the different formula¬ 
tions of the assembly activation event. The three characteristic 
setups we consider are as follows. Figure[6[a) refers to R(D 1 ) 
dependence for the fixed very small £> 2 , whereas in Fig. |6]b) 
is shown the R.{D->) dependence for the fixed very small I)\. 
Finally, Fig. [6]c) illustrates the behavior of R(l) \) when D2 
is fixed at an intermediate value. Comparing Fig. (6ja) and 
Fig. [6jb), it stands out that the formulation 1 on one hand, 
and the formulations 2 and 3 on the other hand, yield sub¬ 
stantially different R dependencies both below and above the 
stochastic bifurcation, cf. Fig. [3}d). Below the bifurcation, 
the respectively smaller R values obtained under formulation 


1 indicate that the corresponding activation process is more 
homogeneous compared to processes conforming to formula¬ 
tions 2 and 3. Furthermore, under formulations 2 and 3, the 
fluctuations around the mean TFP expectedly increase in the 
noise domain that gives rise to stochastic bifurcation. The 
large R values found there indicate strong irregularity of fir¬ 
ing times over an ensemble of different stochastic realizations, 
cf. the peaks in Fig. |7]a) and Fig. [7jb). Sufficiently above 
the bifurcation, the R values decrease for all three formula¬ 
tions. Nevertheless, in the stochastically supercritical state, 
the physical picture regarding the three formulations is in a 
sense reversed compared to the subcritical state. The devi¬ 
ations from the mean TFP with Di are the largest for the 
activation process conforming to formulation 1, whereas R 
values for formulations 2 and 3 are smaller and fairly insensi¬ 
tive to noise increase, at least for sufficient D\. An interesting 
remark is that by adhering to formulation 1, as opposed to for¬ 
mulations 2 and 3, the coefficient of variation never crosses 
1 , the value characteristic for the exponential distribution of 
events. It is well known that the latter distribution is consis¬ 
tent with the Poisson process |36. 381. 

Additional information on R behavior can be gained by 
comparing the TFP distributions of activation events over an 
ensemble of different stochastic realizations for (D 1 , D 2 ) rep¬ 
resentative of the domains below or above the stochastic bifur¬ 
cation (not shown). The main point concerns how likely are 
the deviations biased toward the shorter or longer activation 
times for an ensemble of realizations. Below the bifurcation, 
the ensemble of events under formulations 2 and 3 typically 
splits into fractions with short or very long TFP s. Though 
the former fraction is considerably larger, the long events still 
strongly influence the mean T FP. The partition into two frac¬ 
tions persists under formulation 1, though the difference in 
duration between the two characteristic types of events is con¬ 
siderably reduced. The distribution of TFP s is then found 
to be bimodal, showing two relatively even peaks. These 
points are consistent with the differences displayed by the 
three R(Di,D 2 ) dependencies in Fig[6] in the stochastically 
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FIG. 7. (Color online) Focus on local activation processes, (a) 
and (b) show typical distributions P(ri t i) of individual unit TFP s 
ri : i for a single realization of the assembly activation process at 
(£> 1 , 1 ) 2 ) below and above the stochastic bifurcation, respectively. 
In the subcritical state, one finds two excitation waves, the first being 
induced by noise, whereas the other is due to relaxation of units. In 
the supercritical state, the units are primarily elicited by noise. The 
system parameters are Di = 0., D 2 = 0.00009, c = 0.1, N = 300 
in (a) and Di = 0.0002, D 2 = 0.00025, c = 0.1, N = 300 in (b). 


subcritical state. Nevertheless, above the bifurcation, quali¬ 
tatively similar behavior for all three formulations of activa¬ 
tion event is recovered. As expected, the TFP s follow a uni- 
modal distribution centered around the short events, whereby 
a longer tail becomes visible for formulation 1 if I) \ is in¬ 
creased. The latter accounts for the rise of R with I)\ in Fig. 

Etc). 


B. Distribution of single units’ TFPs 

Having analyzed the properties of collective activation pro¬ 
cess for an ensemble of different stochastic realizations at 
fixed (Di, D 2 ), we now consider the dynamics of local acti¬ 
vation processes for a single realization of an assembly event. 
In particular, we focus on how the TFPs of individual units 
are distributed for noise domains below or above the stochas¬ 


tic bifurcation. As for the underlying mechanism, an interest¬ 
ing point is to examine the contribution of noise-induced local 
activation events relative to the ones elicited by the relaxation 
processes of other units. We find that the typical profiles of the 
distribution of individual TFP s below and above the stochas¬ 
tic bifurcation are quite different, as illustrated in Fig. |7} a) and 
Fig. (TJb), respectively. Note that the qualitatively analogous 
results are obtained if one compares the corresponding distri¬ 
butions of individual T FPs for all units over an ensemble of 
different stochastic realizations (not shown). 

The bimodal distribution in Fig. [71 a) suggests the existence 
of two ’’waves” of local excitation. While the first wave (short 
TFPs) is elicited by noise, the second wave is evoked by re¬ 
laxation of the units that have already emitted the spike (such 
stimuli are conveyed via the interaction terms). This physical 
picture is reminiscent of the one obtained when determining 
the phase response curve of an assembly of oscillators sub¬ 
jected to random perturbations [13^]. Note that the smallest 
TFP s are around 20 a. u., which is about the length of the 
spike excitation loop, i. e. the duration of the typical relax¬ 
ation process. This indicates the existence of activation bar¬ 
rier associated to first pulse emission process of single units. 
Apart from the time scales involved, a remark should also be 
added regarding the sizes of waves, viz. the fractions of units 
participating the waves. The sizes are found to depend on 
D 1 ,D 2 and c, whereby the stronger c expectedly amplifies 
the secondary wave. One would further expect the secondary 
wave to be substantially influenced by the topology of inter¬ 
actions within the network. 

Above the stochastic bifurcation, the local activation events 
no longer maintain the two-wave-like form, see Fig. [7jb). In¬ 
stead, the typical distribution of single unit TFP s becomes 
unimodal, with some differences observed for the cases of in¬ 
termediate vs large noise intensities. In particular, in the latter 
case the peak is more pronounced than in the former, whereby 
the distribution profile itself may be fitted to Lorentzian like- 
form. The very fact that the distribution of local TFP s is 
unimodal in the stochastically supercritical state implies that 
the activation of single units is almost exclusively driven by 
noise. Nonetheless, the point that the TFP s are typically of 
the order of the excitation loop suggests the lack of activation 
barrier in the first pulse emission processes of single units. 


C. Activation mechanisms below and above stochastic 
bifurcation 

While discussing r(D 1 ,D 2 ) behavior in Fig. [3] it has al¬ 
ready been established that below the stochastic bifurcation, 
the assembly activation process under formulations 2 and 3 is 
not primarily affected by the first pulses of individual units, 
but is more likely contributed by the latter spikes. In these 
cases, we have identified gradual coherence build-up between 
individual spikes as the mechanism underlying assembly acti¬ 
vation. Nevertheless, results in Fig. [7jb) suggest that a poten¬ 
tially different scenario of assembly activation is typical above 
the stochastic bifurcation. Therefore, our next goal is to ex¬ 
plicitly demonstrate the distinction between the mechanisms 
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guiding the assembly activation below or above the stochastic 
bifurcation. To this end, we numerically determine the most 
probable trajectories for the first pulse emission process in the 
( X , Y) configuration space and analyze the differences be¬ 
tween the typical paths obtained for ( D \, IJ 2 ) below or above 
the stochastic bifurcation, cf. Fig. [ 8 ] The trajectories pre¬ 
sented conform to formulation 2 of the assembly activation 
event with the fixed threshold value X {) = 0.3. 

The details of the applied numerical method are as follows. 
For the given (Hi, D 2 ), we consider the ensemble of fluctua¬ 
tion paths that start at moment t, from the deterministic fixed 
point (X eq , Y eq ) and reach the terminating boundary defined 
by the threshold X () , cf. the solid line in Fig. [ 8 ] At fj, all the 
individual units are prepared to the same initial conditions, 
i. e. ( x(ti),y(ti )) = (X eq ,Y eq ). The terminating time t f , 
as well as Y(tf) are left unspecified. For the described en¬ 
semble, we consider statistics of the (X(f), Y(t)) position of 
trajectories as a function of time ti < t < tf preceding the 
arrival to Xo = 0.3. Naturally, the recorded trajectories may 
have quite distinct tf times. The proper statistical quantity to 
characterize such an assembly of paths in configuration space 
is the prehistory probability density, first introduced in | 43ll 
and successfully applied many times since [Hi . The former 
can effectively be obtained if the time when each stochastic 
realization terminates is set to t = 0 , such that the behavior 
of the process during the initiation of the pulse is observed by 
looking backward in time. The appropriate prehistory prob¬ 
ability distribution is then defined as H(X,Y,t)dXdY = 
Pr[X(t) e (X,X + dX),Y(t) e (Y,Y + dY)\X(t f ) = 
X 0 ,X(ti) = X eq ,Y (ti) = Y eq ],ti < t < t f ,X < Xq. 
The most probable path for the first pulse emission process is 
determined by collecting the points (X rn (t). Y m (t )) that cor¬ 
respond to the maximum of II (X, Y, t) at the given t. 

In Fig. 0 the most probable activation paths typical for 
the noise domains below (above) the stochastic bifurcation 
are indicated by the empty squares (filled circles). In order 
to explain the differences between the presented paths, it is 
useful to invoke the approximate model ©, which effectively 
assumes that all the units are strongly synchronized. In other 
words, that model refers to scenario where the assembly acts 
as a macroscopic FHN unit. The dashed lines in Fig. [8] de¬ 
note the nullclines of ©. Further, since the latter has the form 
analogous to the dynamics of a single FHN unit, we may 
borrow from our previous paper 121 ] the results obtained for 
such a setup. In particular, we plot the canard-like trajectory 
belonging to the ’’ghost-separatrix” [[ 20 j], cf. the dotted line in 
Fig. [ 8 ] as well as the trajectory generated by the system of 
effective Hamiltonian equations associated to © as discussed 
in our previous paper [21] (see the dash-dotted line in Fig. 
0 . Note that the most probable activation paths for a single 
FHN unit are only sensitive to the external/internal character 
of noise, but not to the particular noise intensities lf 20 l . 21 * 1 . 

For the trajectory corresponding to the noise domain below 
the stochastic bifurcation (empty squares) one finds that the 
approximate model © does not apply. From the latter point, 
one further infers that the units are not coherent during the 
initial part of the first pulse emission process, which is consis¬ 
tent with the already described scenario of gradual build-up 



FIG. 8. (Color online) Most probable activation paths for the 
global variables ( X , Y) under formulation 2 of the activation event. 
The two illustrated profiles are typical for noise domains below 
(empty squares) and above (filled circles) the stochastic bifurcation. 
The shown paths refer to particular noise intensities (£> 1 , 1 ) 2 ) = 
(0,0.00009) and (£> 1 , 1 ) 2 ) = (0,0.0004), respectively. Note that 
the assembly's maximum likelihood trajectories correspond to the 
peak of the histogram obtained for an ensemble of different stochas¬ 
tic realizations as a function of time. The paths for the assembly are 
compared against the one for the approximate model © (’’macro¬ 
scopic FHN unit”). The nullclines (NULL) of this model are 
presented by the dashed lines, whereas the canard-like trajectory 
(CNRD) is denoted by the dotted line. The trajectory generated 
by the effective Hamiltonian system associated with © is given by 
the dash-dotted line. The solid line indicates the threshold Xq = 0.3. 
The results for the assembly refer to system size N = 100 and are 
obtained by averaging over a 1000 different stochastic realizations of 
the first pulse emission process. 


of spike coherence. The mechanism of first pulse emission is 
quite different for the noise domain above the stochastic bifur¬ 
cation. Recall that the stochastic bifurcation underlies transi¬ 
tion to noise-induced oscillations, which are facilitated by the 
(stochastic) synchronization of individual units. With this in 
mind, it is expected to find that the corresponding most prob¬ 
able first pulse emission trajectory (filled circles) conforms 
much better to physical picture behind model ©. In the ini¬ 
tial part of the activation trajectory, there is (stochastic) syn¬ 
chronization between the units of the assembly, such that the 
assembly’s activation path indeed fits well to the path given 
by the ’’macroscopic FHN unit” model ©. Note that the de¬ 
viation from the path corresponding to © is observed in the 
latter part of the assembly activation trajectory. This is a natu¬ 
ral consequence of the fact that the spike times of single units 
in the assembly are synchronized approximately , but are not 
exactly matched. Thus, the effect of noise in the latter part 
of activation trajectory is felt more strongly for an assembly 
than in case of the effectively single FHN unit model ©. As 
expected, the described physical picture on the supercritical 
state breaks down if the noise intensity is too large. 
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FIG. 9. (Color online) Focus on the anti-resonant system-size ef¬ 
fect for r. (a) and (b) illustrate the t(N) dependencies typi¬ 
cal for subcritical noise intensities, whereby (a) conforms to for¬ 
mulation 2 and (b) to formulation 3 of the assembly activation 
event. Both plots are obtained for the same set of parameter val¬ 
ues (D 1 ,D 2 ,c) = (0.0001365, 0.0002255, 0.1), with the threshold 
in (a) fixed to A'o = 0.4. The stochastic averaging at each point has 
been carried out over an ensemble of 1000 different activation paths. 


D. Mean TFP anti-resonance as a system size effect 

In this subsection, we briefly report on a system-size ef¬ 
fect which discriminates between the different formulations 
of the assembly activation event. In particular, the effect is 
observed when applying formulations 2 and 3 which concern 
the dynamics of global variables, but is absent for formula¬ 
tion 1 that refers to local activation events. The system size 
effect consists in the appearance of an anti-resonant peak in 
the dependence of mean TFP with N for the fixed parameter 
set (D i, £> 2 , c). The term anti-resonance is applied in a sense 
that under variation of N, the mean TFP s exhibit a maxi¬ 
mum for the particular size of the assembly. Note that the 
effect is found for (1) \. D->) values both below and above the 
stochastic bifurcation, but is substantially more pronounced 


in the stochastically subcritical state. As an illustration, in 
Figs. |9ja) and [9}b) the t(N) dependencies are shown corre¬ 
sponding to formulations 2 and 3, respectively, whereby the 
noise intensities are typical for the domain below the stochas¬ 
tic bifurcation. Apart from the fact that the formulation 1 does 
not admit similar behavior, it also yields t(N) dependence 
that exhibits only a weak decline with N. As already indi¬ 
cated, under formulations 2 and 3 the anti-resonant effect per¬ 
sists above the stochastic bifurcation. Nevertheless, the t(N) 
maximum in this domain reduces with noise and further shifts 
to smaller N. For a comprehensive understanding of the de¬ 
scribed anti-resonant system size effect, one would have to 
carry out an elaborate analysis on a number of issues indepen¬ 
dent on the present study. The focus should primarily lie with 
the mechanisms contributing the long activation events, in¬ 
cluding the details on how the competition between the noise- 
and relaxation-induced local activation processes is affected 
by the system size. 


VI. SUMMARY 

The present study has been aimed at characterizing the ex¬ 
citable behavior and the noise influenced activation process of 
an assembly of class II excitable units whereby each unit is 
subjected to external and internal noise. The analysis is an 
extension of our previous work m concerning the activation 
processes for a single and two interacting FHN elements. In 
conceptual terms, three points can be regarded as most rel¬ 
evant ones. First, we have explicitly demonstrated that an 
assembly of excitable FHN units may itself exhibit macro¬ 
scopic excitable behavior. As a second point, we have estab¬ 
lished the qualitative analogy between the statistical proper¬ 
ties of the noise-driven activation processes for a single FHN 
unit and an assembly. Finally, it has been shown that depend¬ 
ing on the noise intensities, two qualitatively distinct local 
mechanisms may have prevailing effect on the assembly ac¬ 
tivation process. 

The first point has been achieved by implementing the ap¬ 
propriate MF approach, which ultimately yields the approx¬ 
imate model ©. As in case of a single FHN unit, the 
threshold-like dynamics separating the assembly’s spiking re¬ 
sponses from the small-amplitude excitations has been linked 
to the ghost-separatrix. An additional result is that compared 
to excitable dynamics of a single FHN unit for the analo¬ 
gous parameter set, the dynamics of the assembly exhibits 
larger sensitivity to initial conditions in vicinity of the ghost- 
separatrix. 

As for the second point, establishing the qualitative analogy 
in statistical properties of the activation processes of a single 
unit and an assembly is a non-trivial result that could not be 
expected a priori. To facilitate the comparison between the re¬ 
spective activation processes, we have introduced three alter¬ 
native formulations of the assembly activation event, each em¬ 
phasizing different aspects of collective dynamics. Formula¬ 
tion 1 deliberately makes no mention of global variables, such 
that the effects of coherence of individual units’ spikes are left 
aside. At variance with this, formulations 2 and 3 are stated in 
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terms of global variables, such that the former involves an ex¬ 
plicit scalar threshold for the spiking response, whereas the 
latter incorporates an appropriate terminating boundary set 
analogous to the one we adopted for a single and two inter¬ 
acting FHN units Izfll . 

The analysis on statistical properties of assembly activation 
process has been carried out by determining the r(l)\. I)->) 
and R(Di, Df) dependencies for the three alternative formu¬ 
lations of the activation event. All three formulations yield 
qualitatively analogous results for the mean TFPs , but cer¬ 
tain quantitative differences are manifested in the subcritical 
regime. These differences are due to microscopic mechanisms 
behind the assembly activation, and concern the role of spike 
coherence in the onset of activation event. In particular, the 
activation processes under formulations 2 and 3 are mediated 
by the gradual build-up of spike coherence, whereas in case 
of formulation 1 spike coherence is not relevant. 

We have further found that the statistical properties associ¬ 
ated to macroscopic excitable behavior qualitatively resemble 
those for a single unit [21]. The observed qualitative similar¬ 
ity, both with respect to three formulations of the assembly ac¬ 
tivation event and compared to the case of a single FF[N unit, 
has been tied to stochastic bifurcation that underlies transition 
from the stochastically stable fixed point to continuous oscil¬ 
lations. The stochastic bifurcation has been analyzed within 
the framework of the MF model ([5]i, demonstrating that the 
latter undergoes Hopf bifurcation for the noise intensities that 
qualitatively coincide with those that give rise to stochastic bi¬ 
furcation of the exact system. The stochastic bifurcation has 
been shown to account for the transition between the (D \, !)■>) 
domains admitting large mean TFPs and the plateau region. 

As for our third main point, the analysis of the mechanisms 
behind the local activation processes has revealed substantial 
differences between the stochastically subcritical and the su¬ 
percritical state. In the former case, two excitation waves may 
be discerned: the first is elicited by noise, whereas the second 
wave is due to relaxation of units, viz. is evoked by the in¬ 
teraction terms. For this scenario, one is able to confirm the 
existence of a potential barrier associated to single unit acti¬ 
vation. At variance with this, above the stochastic bifurcation, 
the local activation processes are strongly influenced by noise. 

The respective roles of noise and interaction terms in the 
mechanisms leading to assembly activation have further been 


examined by considering the most probable activation paths 
characteristic for the noise domains below and above the 
stochastic bifurcation. The analysis is based on comparing 
the results for the assembly to the appropriate most probable 
activation path of the approximate model <[6]». The latter model 
assumes strong synchronization between the units, effectively 
describing the assembly as a macroscopic FFIN unit. For the 
subcritical state, the analysis corroborates the gradual build¬ 
up of spike coherence as the leading mechanism behind the 
assembly activation. This mechanism is naturally influenced 
by the interaction terms. Nevertheless, above the stochas¬ 
tic bifurcation, noise facilitates stochastic synchronization be¬ 
tween the units, rendering the interaction terms negligible. In 
this scenario, the assembly activation process may indeed be 
compared to that of the approximate model (|6j. 

In terms of details specific for the assembly excitable be¬ 
havior, an interesting finding concerns a system-size effect, 
which consists in the appearance of an anti-resonant peak in 
the t(N) dependence under fixed (Di, D 2 ). While t(N) dis¬ 
plays a maximum for arbitrary noise intensities, the maximum 
is still substantially more pronounced below than above the 
stochastic bifurcation. 

Apart from providing insights into the assembly activation 
process, the present study has raised a number of novel is¬ 
sues. For example, the future research may include a system¬ 
atic study on the mechanisms behind the reported system-size 
effect, or could focus on the impact of connection topology 
on the activation process. Another relevant point would be 
to consider whether the qualitative aspects of behavior found 
here are paradigmatic, i. e. whether they persist if the assem¬ 
bly is built of type I instead of type II excitable units. 
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